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FOREWORD 


The  work  reported  herein,  covering  the  period  10  April  1972  to 
31  December  1972,  was  carried  out  by  the  Infrared  and  Optics  Division  of 
the  Environmental  Research  Institute  of  Michigan  (formerly  the  Willow 
Run  Laboratories  of  The  University  of  Michigan),  Ann  Arbor,  Michigan. 

The  work  was  performed  under  Contract  DAAD05-72-C-0216  for  the  Army 
Ballistic  Research  Laboratories,  and  was  done  in  three  parts,  each  of 
which  represent  one  volume. 

Toe  three  volumes  are: 

I  -  Polarized  Bidirectional  Reflectance  With 

Lambertian  or  Non-Lambert ian  Diffuse  Component. 

II  -  Polarized  Spectral  Emittance  From  4  to  14  pm. 

III  -  Wavelength  Dependence  of  Polarized  Bidirectional 

Reflectance. 
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TABLE  OF  CONTENTS 


FOREWORD. 


TABLE  OF  CONTENTS. 


LIST  OF  FIGURES. 


LIST  OF  TABLES. 


INTRODUCTION. 


MODEL  DESCRIPTION. 


3 .  MEASUREMENTS . . .  11 

3.1.  Instrument . 11 

3.2.  Spectral  Emittance. . . . 14 

3.3.  Polarized  Emittance  vs.  Angle . . .  19 

4.  MODEL  VALIDATION .  32 

5.  CONCLUSIONS . . . 42 

APPENDIX  A:  DETERMINATION  OF  EXPERIMENTAL  ERROR  IN  MEASUREMENTS 

OF  DEGREE  OF  POLARIZATION .  43 

APPENDIX  B:  DOCUMENTATION  FOR  CLASSICAL  OSCILLATOR  FITTING 

PROGRAM .  45 

APPENDIX  Cl  DOCUMENTATION  FOR  PROGRAM  'EMISPOL' .  .  49 

APPENDIX  D:  PROGRAM  LISTING  -  OSCILLATOR .  51 

APPENDIX  E:  PROGRAM  LISTING  -  EMISPOL .  55 

APPENDIX  F:  PROGRAM  LISTING  AND  INSTRUCTIONS  FOR  IBM 

SCIENTIFIC  SUBROUTINE  -  FMFP .  57 

REFERENCES .  67 

DISTRIBUTION  LIST .  69 


Pnctliit  PW  Wirt 


5 


LIST  OF  FIGURES 


1.  Schematic,  of  FISR  Fore-Optics . . .  12 

2.  Spectral  Directional  Emittance  -  AO-2017 .  16 

3.  Spectral  Directional  Emittance  -  AO-2022 .  17 

4.  Spectral  Directional  Emittance  -  AQ-2023 .  18 

5.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2017  X  =  9.3pm .  20 

6.  Relative  Polarized  Emittance  vs.  Polar  Angle 

aO-2017  A  =  10.6pm .  21 

7.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2022  A  =  5.3pm . . .  22 

8.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2022  A  =  8.3pm . . ,  23 

9.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2022  A  =  9.7pm .  24 

10.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2022  A  =  10.6pm .  25 

11.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2022  A  =  12.0pm . 26 

12.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2023  A  -  5.3pm .  27 

13.  Relative  polarized  Emittance  vs.  Polar  Angle 

AO-2023  A  =  8.3pm .  28 

14.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2023  A  -  9.7pm . . .  29 

15.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2.023  A  -  10.6pm .  30 

16.  Relative  Polarized  Emittance  vs.  Polar  Angle 

A0-2023  A  *>  12.0pm . 31 

17.  Theoretical  Degree  of  Polarization  and  Normal  Emissivity 

vs.  Wavelength  for  O.D.  Paint  (Sample  A01792) .  33 

18.  Model-Measurement  Comparison  Relative  Polarized  Emittance 

vs.  Polar  Angle  A0-2017  A  ■  10.6pm.... .  39 

19.  Model-Measurement  Comparison  Relative  Polarized  Emittance 

vs.  Polar  Angle  AO-2022  A  -  10.6pm .  40 

20.  Model-Measurement  Comparison  Relative  Polarized  Emittance 

vs.  Polar  Angle  AO-2023  A  -  10.6pm...... .  41 

6 


3tPW 


LIST  OF  TABLES 

1.  Polarization  Field  Measurements.. .  34 

2.  Comparison  of  Theoretically  Calculated  and  Experimentally 
Measured  Degree  of  Polarization  for  Emittance  for  X  ■  9.83pm 

(O.D.  Paint  A01792) .  37 

3.  Comparison  of  Theoretically  Calculated  and  Experimentally 

Measured  Degree  of  Polarization  for  O.D.  Paints . .  38 


7 


1.  INTRODUCTION 


Most  materials  on  the  earth's  surface  emit  a  large  portion  of  their 
radiated  energy  in  the  thermal-infrared  wavelength  region.  If  irregularities 
of  the  surface  of  a  natural  target  are  very  large  compared  to  the  wavelength 
of  emitted  radiation,  total  emissivity  measurements  are  the  prime  source 
of  information.  However,  if  the  surface  irregularities  of  a  target  are 
small  compared  to  the  wavelength  of  emitted  radiation,  the  surface  tends 
to  be  more  specular  and  the  electric  vectors  of  this  emitted  radiation 
will  vibrate  in  preferential  directions,  giving  rise  to  an  observable 
emission  polarization.  For  relatively  smooth  surfaces,  therefore,  polar¬ 
ization  measurements  are  an  additional  source  of  information.  Since  the 
surfaces  of  military  targets  are  generally  smoother  and  have  more  ordered 
geometric  patterns  than  naturally  occurring  materials,  polarization  is  a 
potentially  important  parameter  for  target  discrimination. 

The  purpose  of  this  report  is  to  describe  a  phenomenological  model 
for  predicting  emission  polarization  of  paints  using  as  input  only  a 
single,  near-normal  reflectivity  spectrum. 

2.  MODEL  DESCRIPTION 

For  the  case  of  an  emitting  polished  surface,  that  portion  of  the 
emitted  radiance  at  some  angle  d  to  the  target  surface  normal  that  has 
its  electric  vector  vibrating  perpendicular  to  the  plane  of  emission 
(containing  the  surface  norma?  and  the  propagation  vector)  ,  divided  by 
half  the  radiance  emitted  by  a  blackbody  at  the  same  temperature,  is  called 
the  perpendicular  component  F.j^(0)  of  the  emissivity.  [Note:  All  the  emissivity 
components  are  wavelength  dependent,  but  the  \  notation  will  be  suppressed.] 
Likewise,  the  parallel  enissivity  component  E||(6)  defines  that  portion  of 
radiation  emitted  with  electric  vector  vibration  parallel  to  the  plane  of 
emission.  The  total  emissivity  E(0)  is  half  the  sum  of  Ej^(0)  ant*  Ej|(®)» 

The  degree  of  polarization  for  emission  can  be  defined  as: 

PE  -  [E  j  |  (0)  -  Ej^(9)  J/[E  |  j  (9)  +  Ej^(e)  ]  .  (1) 
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The  relationship  between  the  emission  polarization  components  and 
reflection  polarization  components  can  be  found  from  a  consideration  of 
the  reflection  experiment  where  source  (at  0  )  and  observed  (at  0  )  are 
located  on  opposite  sides  of  the  target  surface  normal  at  ar.  angle 

0  =  °i  =  0r*  Let  RJ_  and  R|  |  be  the  ref Activity  components  for  radiation 
with  electric  vectors  vibrating  perpendicular  and  parallel,  respectively, 
to  the  plane  of  incidence  containing  the  incident  and  reflected  rays.  With 
the  constraint  that  9  =  0.  =  0^  implying  that  the  plane  of  emission  coin¬ 

cides  with  the  plane  of  incidence,  Kirchhoff's  law  dictates  (for  equilibrium 
conditions)  that 


E(0)j^=  1  -  RfG)^,  (2a)  E(6)||  =  1  -  R(0 )  j  |  -  (2b) 


The  reflectivity  flux  components  are  functions  only  of  0  =6  -  6 

i  r  ’ 

and  the  complex  index  of  refraction  of  the  target  material  is  given  by 
N  =  n  -  ik.  Expressions  for  R,  and  R. ■  in  terms  of  Ufi,  n,  and  k  are  as 
follows  [1]: 


Rl(6) 


Rj  j  (0) 


[  (a  -  cos  6  > 2  *■  b4'  J 
[ (a  +  cos  0 ) 2  +  b2 ] 


j  ( (a  -  cos  0)2  +  b2] 

[ (a  -  Sine  tan  6 )2  +  b2j 

(  [  (a  +  cos  0 ) 2  +  b2 ] 

X 

2  2 

[ (a  +  sin8  tan  0 )  +  b  ] 

(3a) 


(3b) 


where 


and 


2 

a 


I?  2  9  ^  O  n  n 

2  fn  “  k  “  '’in  '0  +  [4n  k  +  (n  -  k 

•|  {-n2  +  k2  +  sin20  +  [4n2k2  +  (n2  -  k2 


sin2  0  )2)I) 
sin2  0)2F) 


From  Eqs.  (1),  (2),  and  (3),  therefore,  it  is  possible  to  express  the  degree  of 
emission  polarization  P£  as  a  function  of  0,  n,  and  k. 
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For  materials  which  exhibit  reststrahlen  bands  in  the  thennal-IR 

wavelength  region,  n  and  k  are  strongly  wavelength-dependent  in  that 

region.  The  wavelength  dependence  of  P  for  a  given  material  can  be 

r * 

calculated  from  Eqs.  (1),  (2),  and  (3)  if  the  dispersion  curves  of 
n  and  k  are  known.  Once  these  dispersion  curves  have  been  determined 
as  continuous  functions  of  wavelength  for  a  given  paint,  it  is  then 
possible  to  calculate  the  emission  polarization  for  all  6  as  a 
function  of  A.  This  will  be  calculated  for  several  O.D.  paints. 

One  means  for  determining  the  indices  of  refraction  is  the  classical 
oscillator  fitting  method.  The  complex  index  of  refraction  N  -•  n  -  ik 
can  be  calculated  as  a  function  of  frequency  v,  according  to  a  classical 
oscillator  model  of  crystals  [2],  from  the  following  equations: 


2nk 


r  y  .  v  .  v 

yw  v  — 2 — H — 

j  3  3  L(vj  -  v  >  +  Y j 


2  2  2 
V- 


2  ,2 

n  -  k  =■  e,  + 


2  r  v 2  -  v2  -| 


(4) 


where  p^ ,  ,  and  are  the  strength,  width,  and  frequency,  respectively, 

of  the  lattice  oscillator  and  ttj  is  the  high  frequency  dielectric  con¬ 
stant.  If  these  oscillator  parameters  are  known,  n  and  k  are  calculable 
from  Eq.  (4).  If  they  are  unknown,  they  can  be  estimated  by  the  following 
procedure.  Initial  guesses  for  the  oscillator  parameters  are  assumed,  (see 
Appendix  D)  and  a  normal  incidence  Fresnel  reflectivity  curve  R(v)  is  cal¬ 
culated  from  the  equation: 


R(v)  - 

I  (n(v)  +  1]  +  k‘(v)  ) 


(5) 
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A  computer  program,  which  employs  an  IBM  scientific  subroutine,  [see  Appendix 
F]  compares  the  theoretical  reflectivity  curve  calculated  from  Eq.  (5)  with  an 
experimental  normal  incidence  spectral  reflectivity  curve  for  a  smooth  sur¬ 
face  of  the  specimen,  changes  the  oscillator  parameters  in  Eq.  (4), 
calculates  a  new  theoretical  curve,  and  reiterates  until  a  good  fit  is 
made  between  theoretical  and  experimental  curves.  With  the  oscillator 
parameters  of  the  best-fitting  theoretical  curve,  the  indices  can  then 
be  calculated  from  Eq.  (4).  (See  Appendix  D  for  use  of  Oscillator  Program.) 

3.  MEASUREMENTS 

The  measurements  in  this  program  consisted  of  two  phases.  The  first 
phase  consisted  of  a  spectral  scan  of  the  emittanca  from  each  sample  with 
a  normal  angle  of  observation.  The  second  phase  consisted  of  angular 
scans  of  the  emittance  with  receiver  polarizer  set  alternately  perpendicular 
and  parallel  to  the  receiver  incidence  plane.  In  this  latter  case,  each 
scan  took  place  at  a  particular  wavelength. 

In  the  following  paragraphs  we  provide  a  brief  description  of  the 
instrument  which  was  used  for  these  measurements,  followed  by  data  which 
resulted . 

3.1  The  Instrument 

The  field  infrared  spectroradiometer  (FISR)  which  was  utilized  is  a 
double-beam  Instrument  with  output  proportional  to  the  radiance  difference 
between  the  objects  filling  the  fields  of  the  two  Herschelian  reflecting 
telescopes.  The  measurement  of  the  difference  is  accomplished  by  a  reflective 
chopper  which  alternately  samples  the  radiation  from  the  two  telescopes. 

After  being  chopped,  the  radiation  fro.i  the  two  telesceopes  is  imaged  at  a 
common  focal  plane  at  which  a  field  stop  slit  is  placed,  limiting  the  FOV 
to  a  0.6‘  x  2.5*  field.  The  fore-optics  just  described  are  illustrated  in  Figure 
1.  The  chopped  energy  passing  through  the  field  stop  slit  is  refocused, 
by  a  90s  off-axis  ellipsoid  mirror,  onto  an  entrance  aperture  which  lies 
immediately  in  front  of  the  circular  variable-filter  (CVF)  monochromator. 
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Figure  1.  SCHEMATIC  OF  FI°R  FOItE-OFTlCS 


,  12 


The  CVF  then  transmits  a  narrow  spectral  band  of  this  energy  to  a  45°  off- 
axis  ellipsoid  mirror,  which  refocuses  it  onto  the  detector. 

The  generated  ac  detector  signal  is  amplified  and  then  synchronously 
demodulated.  The  synchronous  demodulator  transforms  those  components  of 
the  signal  which  have  the  particular  frequency  of  the  chopper  modulation 
into  a  proportional  dc  signal,  while  rejecting  those  components  with 
different  frequency  (i.e.,  noise).  Thus,  at  the  output  a  dc  signal  is 
presented  which  has  a  value  at  any  point  in  time  proportional  to  the  energy 
difference  between  the  two  optical  beams  at  a  particular  wavelength  deter¬ 
mined  by  the  CVF  position.  Mechanical  rotation  of  the  CVF  provides  a 
spectral  scan.  The  output  voltage  is  recorded  by  using  a  digital  voltmeter 
as  an  analog-to-digital  converter  and  printing  on  a  digital  paper-tape 
printer,  and  the  wavelength  is  registered  simultaneously  by  printing  the 
digital  output  of  a  shaft  encoder  mechanically  linked  to  the  CVF. 

The  spectrometer  portion  of  the  system  utilized  a  circular  variable 
filter  for  spectral  dispersion  in  the  range  from  2.5  to  14.5pm.  The  CVF 
has  the  advantage  of  providing  more  efficient  energy  transmission  than  a 
prism  disperser,  allowing  faster  spectral  scan  rates  or  narrower  spectral 
resolution  for  a  given  energy  input. 

In  order  to  achieve  high  system  sensitivity,  an  indium  antimonide  and 
a  mercury  doped  germanium  detector  are  used  to  cover  the  wavelength  regions 
2.5  to  5.5pm  and  4.5  to  14.5pm,  respectively.  These  detectors  are  readily 
interchangeable,  so  one  may  be  replaced  by  the  other  during  the  course  of 
scanning  from  2.5  to  14.5pm  without  stopping  the  instrument. 

Two  temperature-controlled  blackbody  calibration  sources  are  provided 
for  insertion  at  the  entrance  optics  in  either  or  both  of  the  two  channels. 
In  the  normal  mode  of  operation,  a  blackbody  remains  in  one  channel,  and 
the  other  channel,  having  previously  been  calibrated,  views  the  target  of 
interest.  In  this  way  the  system  can  be  used  to  measure  absolute  radiance. 
Alternatively,  if  the  target  temperature  is  known,  or  if  a  black  surface 
at  the  same  temperature  is  available,  the  spectral  emittance  of  the  surfae 
can  be  determined.  These  blackbodles  are  only  suitable  for  the  spectral 
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range  4.5  to  14.5pm.  Calibration  sources  below  4.5pm  have  not  as  yet 
been  implemented.  However,  for  these  reflective  wavelengths  (less  than 
4.5pm),  a  calibrated  reflectance  panel  may  be  inserted  in  the  reference 
beam  and  the  instrument  then  used  to  make  relative  reflectance  measure¬ 
ments. 

Supplementary  instrumentation  has  been  developed  to  permit  the  basic 
FISR  equipment  to  be  used  for  obtaining  directional  polarized  emissivity 
data.  Two  similar  samples  maintained  at  different  known  temperatures  are 
separately  viewed  through  each  port  of  the  FISR,  and  the  difference  in 
radiance  values  is  measured.  From  this  measured  difference,  with  the 
available  temperature  information,  the  emissivity  can  be  calculated.  The 
sample  temperatures  are  maintained  by  thermally  controlled  water  baths. 

The  sample  holder  is  designed  to  be  operated  either  at  a  fixed  viewing 
angle  or  by  revolving  and  permitting  a  scan  of  the  sample  surface. 

Complete  details  of  the  FISR  can  be  found  in  reference  4. 

3.2  Spectral  Emlfctance 

Figures  2,  3  and  4  show  the  spectral  emittance  data  in  plotted  form. 

In  each  case  there  are  many  scans,  some  of  which  cover  different  spectral 
regions  than  others.  Each  scan  number  has  its  own  symbol  so  it  can  be 
roughly  identified  on  the  plot.  Shown  also  are  curves  representing  both 
accuracy  and  precision  of  the  measurements.  Accuracy  is  determined  with 
respect  to  both  systematic  and  random  errors.  Precision  is  determined  with 
respect  to  noise  in  the  recorded  voltage,  (Reference  4  provides  a  complete 
discussion  of  determination  of  accuracy  and  precision  in  this  context.) 

Note  that  for  all  three  samples,  there  is  an  apparent  dip  in  the  spectrum 
between  9.5pm  and  10pm.  This  dip  is  interpreted  to  be  a  reststrahlan  band. 
Structure  of  this  sort  arises  from  the  wavelength  dependence  of  the  indices 
of  refraction.  (In  section  1,  we  have  described  how  one  calculates  the 
indices  of  refraction  by  use  of  a  classical  oscillator  fitting  method.) 
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Note  also  that  in  Figure  3  there  is  an  apparent  discrepancy  in  that 
scans  numbers  3  and  4  are  translated  upward  with  respect  to  scans  2  and  23. 

This  discrepancy  appears  to  be  due  to  a  combination  of  calibration  error 
and  difficulty  of  maintaining  accurate  measurement  of  environment  temperature 
during  the  measurements.  The  problem  is  still  being  analysed.  i„  the  validation 
(see  Section  4)  scans  3  and  4  were  used.  In  view  of  the  apparent  difficulty 
it  is  not  clear  that  these  were  the  most  appropriate. 


mm 
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FIGURE  2 


-ENGTHQW 

FIGURE  4 


3.3  Polarized  Emittanee  vs .  Angle 

Figures  5  through  16  show  the  variation  for  each  sample  at  a  variety 

of  wavelengths  of  relative  polarized  emittanee.  In  each  case  the  curve 

E||(e> 

which  represents  parallel  polarization  is  a  plot  of  ej  |  (©)  -  g  and 

e,  <e) 

that  which  represents  perpendicular  polarization  is  ej^(8)  ”  £('())'• 

Therefore,  for  each  angle,  degree  of  polarization  is  extracted  from 
the  data  by  equation  (1),  where  substituting  e | | (6>  and  ei(6)  for  Ej^(0) 
and  E.i(6)  leaves  P_  unchanged,  since  Eii(O)  *  Ei(0)  »  E(0).  [See  Appendix 


A] .  Comparison  of  model  predictions  with  selected  measurement  data  from 
this  group  is  described  in  Section  4. 
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FIGURE  7.  Relative  Polarized  Emittance  vs.  Polar  Angle 
A0-2022  a  «*  5.3um 


FIGURE  11.  Relative  Polarized  Emictance  vs.  Polar  Angle 
A0-2022  A  =  12.0Mm 


FIGURE  12.  Relative  Polarized  Etnictance  vs.  Polar  Angle 
AO-2023  A  »=  S.Jmn 
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4.  MODEL  VALIDATION 


Figure  17  shows  results  of  the  emittance  model  for  an  0.  D. 
paint,  calculated  last  year  as  part  of  the  Target  Signature  Analysis 
Program  of  ERIM,  for  the  U.  S.  Air  Force.  In  this  figure  the  degree  of 
polarization  of  the  emittance  of  an  0.  D.  paint  sample  was  calculated 
as  a  function  of  wavelength  for  the  case  of  no  reflective  input  other 
than  clear  sky  radiance.  The  target  was  assumed  to  be  at  a  temperature 
T  *  300°K  and  the  angle  of  observation  was  chosen  to  be  6^  ■  70° .  In  this 
theoretical  treatment,  the  paint  was  assumed  to  be  a  perfectly  specular 
surface,  with  neglible  volume  and  multi-surface  reflection.  A  normal- 
incidence  reflectivity  curve  was  fit  with  the  classical  oscillator  model 
discussed  above  to  find  values  for  the  complex  index  of  refraction 
N  ■  n  -  ik  as  a  function  of  wavelength.  These  indices  were  then  fed 
into  Fresnel  equations  for  R|  j  and  Rj^  (Eq.  3)  from  which  E j  |  (0)  and  Ej^(b) 
were  calculated  from  Eq.  (2).  Emission  polarization  was  then  calculated 
from  Eq.  (1). 

The  degree  of  polarization  of  most  naturally-occurring  materials  is 
low  in  the  thermal  IR,  compared  to  the  polarization  of  painted  surfaces. 

The  only  known  measurements  of  emission  polarization  from  natural  materials 
are  given  in  Table  I,  which  shows  the  degree  of  polarization  for  a  few 
natural  targets  measured  broadband  (approximately  8  to  14pm)  under  clear, 
humid  sky  conditions  [3].  The  highest  degree  of  polarization  is  approximately 
21  for  these  materials,  as  compared  with  approximately  9X  for  an  O.D.  paint 
sample  measured  at  the  same  angle  (6  ■  75°)  under  the  same  conditions. 

Had  the  measurements  been  in  a  2pm-wide  band  centered  near  9.8pm,  Figure  17 
predicts  that  for  the  O.D.  paint  the  degree  of  polarization  would  have 
been  appreciably  higher.  The  increase  in  the  degree  of  polarization  in 
the  reststrahlen  band  for  target  materials  may  be  useful  in  the  discrimina¬ 
tion  of  such  targets  from  backgrounds  like  water,  which  also  can  have  a 
high  degree  of  polarization. 
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FIGURE  17 


Name 


Packed  Beach 
Sand 

Packed  Beach 
Sand 

Wet  Sand 

Granite  Rock 
(Irregular 
snaps,  nat¬ 
ural  surface) 


TABLE  I 

POLARIZATION  FIELD  MEASUREMENTS 

(August  5,  1970,  Clear 
Sky  Conditions,  Willow 
Run  Airport) 


•  detected  radiance  from  target  with  electric  vector  vibrating 
parallel  to  plana  of  emission 


Lx  ■  detected  radiance  from  target  with  electric  vector  vibrating 
perpendicular  to  plane  of  amission 


For  the  natural  materials  in  Table  I,  the  degree  of  polarization 
generally  increased  with  increasing  observation  angle  for  9  >  50°,  but 
the  measured  polarization  was  so  slight  for  6  <  50°  that  it  fell  below 
noise  levels  of  the  radiometer  system.  As  for  angular  dependence  of 
polarization  from  0.  0.  paint.  Table  II  gives  measured  and  calculated 
values  for  the  degree  of  polarization  of  the  emlttance  for  6  from  10°  to 
70°  for  a  wavelength  of  9.83pm,  the  center  of  the  reststrahlen  band. 

For  6  >  40°  the  emission  polarization  of  paints  can  be  quite  large  com¬ 
pared  with  that  of  natural  target  materials. 

The  good  agreement  between  calculated  results  and  the  measured  values 
shown  in  Table  II  were  encouraging  enough  to  justify  further  validation 
with  respect  to  other  O.D.  paint  samples.  Under  this  BRL  contract,  three 
O.D.  painted  surfaces  were  studied  using  an  ERIM  spectrometer  which  enabled 
us  to  measure  spectral  radiance  distributions  as  well  as  polarization  de¬ 
pendent  radiance  as  a  function  of  observation  angle.  The  samples  are 
labeled  2017,  2022,  and  2023. 

Compared  with  the  earlier  sample  (1792),  these  three  O.D.  paints  ex¬ 
hibited  less  pronounced  reststrahlen  bands;  hence,  their  degree  of  polar¬ 
ization  is  neither  as  large  or  as  wavelength  dependent.  Table  III  shows 
the  measured  and  calculated  (from  the  above  model)  degree  of  polarization 
for  the  three  O.D.  paints  2017,  2022,  and  2023  at  angular  increments  of 
10*  for  wavelengths  of  10.6pm,  12.0pm,  9.7pm,  8.3pm,  and  5. 3pm.  No 
calculations  were  made  for  the  5.3pm  wavelength  primarily  because  the  data 
were  very  noisy  below  7pm. 

As  can  be  observed  from  the  "Model  Result"  sod  "Measured"  columns,  the 
model  always  predicts  a  higher  degree  of  polarization  than  is  measured, 
owing  to  the  assumptions  in  the  model  that  the  paint  surfaces  are  perfectly 
smooth  and  that  the  paint  layer  is  optically  thick.  The  middle  column 
under  each  paint  sample  shows  the  product  of  multiplying  the  sample  result 
times  a  constant  factor  of  0.75.  The  physical  reason  for  choosing  a  constant 
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correction  factor  is  that  the  surface  roughness  should  affect  radiation 
from  all  wavelengths  and  angles  about  the  same,  if  the  roughness  is  on 
a  much  different  scale  of  magnitude  than  the  wavelength  variation  (5pm  - 
12pm)  used  for  the  calculations  of  this  model.  The  magnitude  of  the 
constant  factor,  however,  was  chosen  solely  on  an  empirical  basis.  As 
shown  in  Table  III,  the  product  of  the  model  result  and  the  0.75  factor 
is  within  experimental  error  of  the  polarisation  measurements  for  6  of  6 
measurements  on  0.  D.  paint  2017,  for  22  of  28  measurements  on  O.D.  paint 
2022,  and  for  15  of  27  measurements  on  O.D.  paint  2023.  The  calculation 
of  experimental  error  is  explained  in  Appendix  A.  Overall,  these  cal¬ 
culated  results  for  three  paint  samples  and  4  wavelengths  fell  within 
experimental  error  for  43  of  61  measurements,  or  70%  of  the  time.  On 
this  basis,  the  model  (including  the  0.75  constant  multiplicative  factor) 
has  been  verified  for  degree  of  polarization  calculations. 


Figures  18,  19,  and  20  are  plots  of  relative  polarized  emlttance  versus  > 

polar  angle  at  A  -  10.6pm  for  O.D.  paints  2017,  2022,  and  2023  respectively.  j 

They  show  how  the  parallel  and  perpendicular  components  of  emlttance  cal-  j 

■J 

culated  by  the  model  compared  with  experimental  measurements.  No  multi-  j 


plicative  or  other  factors  have  been  employed  to  alter  the  oscillator  model 
results  in  these  figures.  The  relative  polarized  emlttances  ere 


E  I  (©) 

.  Figure  18  shows  a  good  correlation  bet¬ 


ween  theoretical  and  observed  relative  polarized  emlttances  for  paint  2017, 
where  4  unpolarized,  near-normal  spectral  emlttance  curves  were  averaged 
to  produce  the  input  for  the  oscillator  model.  Figure  19  shows  what  can 
happen  when  the  unpolarized  spectral  emlttance  curves  have  measurement 
discrepancies.  An  average  of  the  best  two  curves  was  used  as  oscillator 
model  input  to  produce  the  dashed  line.  The  two  worst  curves  were  used  to 
produce  the  results  shown  by  the  rectangles,  and  the  circles  show  what 
happens  when  all  four  curves  are  averaged  prior  to  application  of  the 
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TABLE  II 


COMPARISON  OF  THEORETICALLY  CALCULATED  AND 
EXPERIMENTALLY  MEASURED  DEGREE  OF  POLARIZATION 
FOR  EMITTANCE  FOR  X-  9.83pm 
(O.D.  Paint  A01792) 


e 

Calculated 

Degree  of  Polarization 

Measured 

Degree  of  Polarization 

10* 

0.58% 

0.25% 

20* 

2.371 

1.06% 

30® 

5.45% 

3.47% 

!  40® 

10.03% 

7.50% 

50® 

16.35% 

12.56% 

60® 

24.71% 

20.25% 

70® 

35.44% 

32.46% 

■***»  •^ts*ei6aBl(«-.^f!ft*T.* 


TABLE  III 

COMPARISON  OF  THEORETICALLY  CALCULATED  AND  EXPERIMENTALLY 
MEASURED  DEGREE  OF  POLARIZATION  FOR  O.D.  PAINTS 


O.D. 

Paint  AO 2017 

0 

.D.  Paint  AO 20 22 

a 

D.  Paint  AO 20 23 

Calculated 

Calculated 

Calculated 

e 

(X) 

Measured 

(X) 

(X) 

Measured 

Modal 

Multiplied 

(X) 

Model  2 

(ultiplied 

1 

lodel 

Multiplied 

Rasult 

by  0.75 

Result 

by  0.75 

V*/ 

Lesult 

by  0.75 

V*/ 

X  * 

LO* 

0.3 

0.2 

0.341.9 

0.3 

0.2 

0.341.5 

0.6 

0.5 

0.9+1. 8 

10.64m 

20* 

1.4 

1.0 

1.941.9 

1.3 

1.0 

1.041.5 

2.3 

1.7 

I.44I.8 

JO* 

3.3 

2.5 

2.841.9 

2.1 

2.3 

2.241.5 

5.1 

3.8 

2.7±1.8 

LO* 

6.1 

4.6 

4.541.9 

5.7 

4.3 

4.341.5 

9.0 

6.8 

4. 1±1. 8 

50* 

10.2 

7.7 

7.741.9 

9.5 

7.1 

7.041.5 

13.7 

10.3 

7.341.8 

SO* 

15.9 

11.9 

2.342.0 

14.6 

11.0 

11.441.6 

19.2 

14.4 

11.341.9 

70* 

23.7 

17.8 

- 

21.4 

16.1 

16.741.7 

25.5 

19.1 

17.0*2.1 

X  - 

LO* 

0.3 

0.2 

<0.341.9 

0.4 

0.3 

<0.542.0 

12. On® 

20* 

1.1 

0.8 

0.741.9 

1.8 

1.4 

<0.542.0 

JO* 

2.6 

2.0 

2.641.9 

4.1 

3.1 

1.642.0 

10* 

5.0 

3.8 

4.541.9 

7.5 

5.6 

3.342.0 

50* 

8.4 

6.3 

7.541.9 

12.0 

9.0 

6.642.0 

SO* 

13.3 

10.0 

12.042.0 

17.6 

13.2 

10.2*2.1 

70* 

20.0 

15.0 

18.042.2 

24.7 

18.5 

15.7*2.3 

X  - 

LO* 

0.4 

0.3 

<0.342.0 

0.7 

0.5 

<0.542.0 

9.7ua 

20* 

1.7 

1.3 

1.342.0 

2.7 

2.0 

<0.542.0 

JO* 

3.8 

2.9 

2.942.0 

5.9 

4.4 

2. 4+2.0 

L0T 

6.9 

5.2 

5.442.0 

10.0 

7.5 

4.642.0 

SO* 

11.0 

8.3 

9.742.0 

14.7 

11.0 

8.0±2.0 

so* 

16.2 

12.2 

14.642.1 

19.8 

14.9 

12.5*2.1 

70* 

22.6 

17.2 

21.742.4 

25.3 

19.0 

L8.6*2.4 

X  - 

LO* 

0.4 

0.3 

<0.341.8 

0.8 

0.6 

1. 0*1.9 

8. 3|i« 

20* 

1.8 

1.4 

0.741.8 

3.3 

2.5 

1.641.9 

JO* 

4.2 

3.2 

2.041.8 

7.1 

5.3 

3.5±1.9 

AO* 

7.4 

5.6 

3.541.8 

11.5 

8.6 

5.641.9 

50* 

11.3 

8.5 

6.341.8 

15.9 

11.9 

9.841.9 

60* 

L5.4 

11.6 

9.541.9 

20.0 

15.0 

.5. 6*2.0 

70* 

L9.6 

14.7 

L5. 042.1 

23.8 

17.9 

- 

X  • 

10* 

0.2 

- 

S.3u® 

20* 

0.8 

0.3 

30* 

1.8 

1.5 

40* 

3.2 

2.9 

50* 

5.5 

5.1 

sa 

9.0 

* 

8.7 

13.9 

FIGURE  18.  Model'-Meaauraaent  Coaparlaon-Relatlve  Polarized 
Eailttance  va.  Polar  Angle  A0-2017  X  *  10.6ua 
•**■■■  Experimental 
..........  lhttorteic#1 


FIGURE  15.  Model-Measurement  Comparison-Relative 
Polarized  Emittance  vs.  Polar  Angle 

AO-2022  A  -  10.6ym 

■*■■■  Experimental 

........ m  Theoretical  using  best  2  scans  of  spectral 

emittance  data. 

G3  *  Theoretical  using  worst  2  scans  of  spectral 
emittance  data. 

O  »  Theoretical  using  all  4  scans  of  spectral 
emittance  data. 


oscillator  model.  This  points  out  the  need  for  better  spectral  emittance 
measurements  as  Inputs  to  the  model.  Figure  20  shows  the  poorer  results 
with  paint  2023.  The  oscillator  model  was  applied  twice  to  the  average 
of  4  spectral  emittance  curves,  with  almost  the  same  results.  The  poorer 
agreement  in  Figure  20  is  obvious  also  in  the  smaller  than  expected 
calculated  degree  of  polarization  results  of  Table  111  for  paint  2023. 

5.  CONCLUSIONS 

The  oscillator  model  for  calculating  polarised  components  of  spectral 
emittance  and  degree  of  polarization  in  the  thermal  infrared  region  for 
smooth-surfaced  targets  has  produced  theoretical  values  that  are  within 
experimental  measurement  accuracy  701  of  the  time.  The  accuracy  of  the 
model  could  be  improved  with  better  measurements  of  unpolarlzed,  near¬ 
normal  spectral  emittance  values.  Such  measurements  could  be  improved  by 
utilizing  a  parabolic  ref lactometer  in  the  reflection  mode  instead  of  FISR,  in 
the  emission  mode  because  of  the  greater  absolute  accuracy  afforded  by  the 
former.  FISR  measurements  would  still  probably  be  preferred  for  the  re¬ 
lative  polarized  emittance  versus  polar  angle  measurements,  because  of  the 
induced  polarization  inherent  to  the  parabolic  ref lactometer.  The  model 
itself  could  probably  be  improved  by  accounting  for  thln-film  transmission 
of  paint  on  metal,  but  this  might  not  be  required  after  better  spectral 
emittance  measurements  are  available. 

Although  the  O.D.  paints  used  in  this  investigation  are  similar  in 
spectral  emittance  features,  the  oscillator  model  should  work  at  least  as 
well  for  other  smooth-surfaced  materials  which  exhibit 
marked  dispersion  (high  X-dependence  of  the  spectral  emittance)  in  the 
thermal  IR  region,  such  as  teflon,  mylar,  and  possibly  phenolics. 

The  most  important  feature  of  this  model  is  that  the  degree  of  polar¬ 
isation  can  be  calculated  for  smooth-surfaced  materials  In  the  thermal  IR 
region  for  any  X  (from  about  7  to  15pm)  and  any  polar  angle,  with  only  a 
measurement  of  spectral  emittance  from  7  -  15pm  at  near-normal  Incidence 
required.  This  will  save  considerable  time  and  money  on  experimental 


measurements . 
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APPENDIX  A 

DETERMINATION  OF  EXPERIMENTAL  ERROR  IN 
MEASUREMENT  OF  DEGREE  OF  POLARIZATION 


The  equation  for  degree  of  polarization  can  be  written  aa: 

E  j  j<0)  -  E_L<6)  e||(e)-e|(e) 

PE(0)  "  E  j  |  (0)  +  Ej^(8)  "  e|  j  (0)  +  ej^(0)  (A-1) 

E(0)|| 

where  e| | (6)  -  E^0j‘  -  relative  polarized  emittance,  parallel  component 

E(6)_L 

e. (0)  a  -  -  relative  polarized  emittance,  perpendicular  component 

1  E(0) 


All  of  the  above  and  following  parameters  are  wavelength  and  angle  dependent, 
but  the  \  and  0  notation  will  be  suppressed  here. 

The  error  in  P£  due  to  errors  in  measured  values  of  e  1 1  and  respectively, 
can  be  represented  by: 


1 


'll  +  el 


'll -'1  1 

('ll  +'l)2  J 

'll  ~'l  1 

"ir*i,2J 


2ei 


( <e | |  +  ei ) 


-2e 


+  ei> 


(A-2) 


(A-3) 


The  total  error  in  PA  can  be  described  as: 


AP- 


(AP  |  |  >  2  ■*"  (APj^)  2  - 


(®| |  +  «j) 


/  2  2 
/(ej^Aejj)  +  ( e j  j Aej^) 


(A™4) 


Let  o  be  the  precision  error  in  the  measurement  of  unpolarised  spectral 

r 

emittance.  It  will  be  assumed  that  Aej |  and  Aej^  are  equal  and  that  the 
addition  of  a  polarizer  into  the  optical  system  increases  the  precision 
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error  by  a  factor  of  /2,  such  that 

Ae  1 1  ■  Aej^  »  (A-5) 

When  equation  (A-5)  is  substituted  into  equation  (A-6) ,  the  resultant  error 

I 


in  Pg  from  error  measurements  in  e  i  i  and  e  |  is : 


2  o 


APE“ 


(•|  |  +  mj) 


y  /2(e2j  |  +  ej2)  -  2.828  0j 


/2  .  2 

/e  M  +  e 


I 


(ej  |  +  ej) 


(A-6) 


TheAPgis  multiplied  by  100  in  Table  Ill  to  get  X  APg.  It  can  be  observed 
from  equation  (A-6)  that  APg  is  not  a  strong  function  of  the  values  of  ej | 
and  ej~,  and  thatAP^^  2.828ap  for  small  9,  where  e|  |  and  e^  are  near  1.0 
in  value. 
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APPENDIX  B 


DOCUMENTATION  FOR  CLASSICAL 
OSCILLATOR  FITTING  PROGRAM 

The  oscillator  fitting  program  is  designed  to  compute  oscillator 
parameters  from  which  effective  indices  of  refraction  can  be  calculated, 
according  to  Eq.  (1)  of  the  text.  Initial  estimates  of  the  number  of 
oscillators  responsible  for  spectral  features  of  the  specimen  and 
the  strengths  (S),  frequency  positions  (NU  ■  iq^qqq,  for  *  in  pm),  and 
widths  (GAMMA)  of  each  oscillator  must  be  made.  The  dielectric  constant 
at  infinite  frequency  (E)  must  also  be  made  known.  [Note:  Here  E  *»  Ew.] 

To  begin  with  one  assumes  at  least  one  oscillator  is  located  at 
each  reflectance  maximum  in  the  Infrared  reflectance  spectrum.  This 
determines  the  initial  number  cf  oscillators  and  their  positions  (NU) . 

The  initial  values  of  E,  S,  and  GAMMA  can  be  taken  to  be  equal  to  the 
values  determined  for  the  O.D.  paint  2017  in  Table  B-l.  Table  B~1  gives 
the  result  of  the  oscillator  fitting  program  for  the  three  O.D.  paints 
listed  in  Table  1  of  the  text. 

The  program  will  refine  the  position,  strength,  and  shape  of  the 
given  oscillators,  but  it  will  not  add  or  delete  oscillators.  However, 
if  an  oscillator  is  grossly  superfluous,  the  strength  will  be  decreased, 
and/or  the  width  will  be  increased,  to  such  an  extent  as  to  make  its 
contribution  negligible. 

The  program  returns  an  error  parameter  before  supplying  the  new 
oscillator  parameters.  If  the  error  parameter  is  0,  convergence  (under 
the  input  EPS  and  EST?  has  taken  place,  otherwise  it  has  not.  Whether 
or  not  convergence  has  taken  place,  it  is  often  advisable  to  look  at  the 
oscillator  parameters  and  the  emittance  values  calculated  from  them  to 
decide  if  additional  oscillators  should  be  added  and/or  subtracted,  or 
whether  the  liming  conditions  (ESP,  EST)  should  be  changed.  If  the 
theoretical  and  (experimental  reflectance  curves  do  not  quantitatively 
match,  i.e.  do  not  peak  and  trough  near  the  same  wavelengths,  addition  or 
subtraction  of  an  oscillator  is  called  for.  In  general,  two  things  are 

*See  last  page  of  this  Appendix  for  definitions. 
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necessary  to  judge  a  good  fit:  a  sun  of  the  differences  squared  (between 

_2 

theoretical  and  experimental  spectral  reflectance)  less  than  about  10 

and  a  qualitative  match  between  peaks  and  valleys  of  experimental  and 

theoretical  spectral  reflectance  curves.  For  low-contrast  materials 

(small  reststrahlen  bands),  a  good  quantitative  fit  is  usually  not  a  sufficient 

criterion  for  successful  fitting,  as  seen  by  the  results  for  sample  2023 

in  Table  B-l  and  Table  1  of  the  text.  (See  discussion  in  text  of  Appendix  C.) 
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TABLE  B-l 

EFFECTIVE  OSCILLATOR  PARAMETERS  FOR  POINTS  RESULTING 
FR OM  CLASSICAL  OSCILLATOR  FITTING  PROGRAM 


0.  D.  Paint  Sample 


Oscillator  Parameters 


E 

S 

NU 

GAMMA 

AO 2017 

1.040 

0.5421 

0.1053 

993.7 

1373.9 

0.3234 

0.4264 

A02022 

1.051 

0.0866 

0.5716 

787.8 

988.4 

0.1190 

0.4168 

AO 202 3 

.726  ! 

1.076 

828.0 

.7384 

j 

.2456 

1022.6 

4.683 

mm 
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DATA  DECK  FOR  THE  OSCILLATOR  PROGRAM 


CARD  »  FORMAT 

1  20A4 


2  15,  E15.5, 

and  F10.0 


3 


15 


4 


3F10.0 


VARIABLES (a) 
A 


LIMIT, EPS ,EST 


N 


S.NU, GAMMA 


DESCRIPTION 

Any  title  up  to  80 
characters  in  length 
(Including  blanks). 

LIMIT  »  number  of  interations 
desired  (30  is  a  good  choice); 
EPS  ■  test  value  for  the 
expected  absolute  error 
(see  the  write-up  on  FMFP; 
despite  what  it  says  0.05 
has  worked  reasonably  well) ; 
EST  “  estimated  value  of  F(x) 
(also  see  Appendix  F). 

Three  times  the  number  of 
cards  of  type  4,  plus  1  (i.e., 
total  number  of  parameters 
needed  to  describe  all  the 
oscillators  with  card  types 

4  and  5) . 

5  »  strength  of  the 
oscillator;  NU  ■  location 
of  the  osclllaotr  (in  wave 
numbers) ;  GAMMA  -  width 
of  oscillator. 


5  F10.0  E 

6  15  M 

7  2F10.0  WIN(I) 


Dielectric  constant  at  in¬ 
finite  frequency 

Number  of  points  used  to 
describe  the  curve  (i.e., 
no.  of  cards  of  type  7). 

RIN(I) 

WIN(I)  •  wave  number  for  the 
data  point;  RIN(I)  •  percent 
emittance  at  WIN(I).  Use  as 
many  of  these  as  necessary  to 
throughly  delineate  the  curve 
(at  least  45-50). 
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APPENDIX  C 


DOCUMENTATION  FOR  PROGRAM  'EMISPOL*. 

EMISPQL  is  a  straightforward  program,  in  which  the  oscillator  para¬ 
meters  from  the  previous  program  are  input  parameters.  Output  is  degree  of 
polarization  (PERPOL)  and  e | j (8)  and  a^(0) ,  relative  emittance  components. 

Where  reststrahlen  bands  are  small,  it  is  important  to  note  that  the 
index  of  refraction  as  determined  by  the  oscillator  program  is  not  reliable. 
If  band  structure  is  not  pronounced  it  is  difficult  to  determine  suitable 
oscillator  parameters.  On  the  other  hand,  the  emissivity  is  then  not  very 
spectrally  dependent,  which  means  that  it  is  reasonable  to  regard  the 
refractive  index  as  constant  over  the  same  wavelength  region. 

Under  these  circumstances  it  is  valid  to  by-pass  the  oscillator  pro¬ 
gram  and  feed  a  constant  value  for  refractive  index  into  EMISPOL.  Doing 
so  requires  a  slight  modification  of  EMISPOL.  Inanedlately  after  statement 
600  (see  EMISPOL  listing  in  Appendix  E) ,  the  program  calls  in  a  subroutine 
which  calculates  n  and  k  given  the  final  oscillator  parameters.  This 
subroutine  (and  hence  the  oscillator  program)  can  be  by-passed  by  sub¬ 
stituting  a  statement  which  supplies  values  for  n  and  k.  Values  can  be 
determined  from  the  Brewster  angle  as  described  in  Volume  I  of  this  report. 
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DATA  DECK  FOR  EMISPOL 


CARD  # 

FORMAT 

NAME 

DESCRIPTION  AND/OR  COMMENTS 

1 

10A4 

ITITLE 

Any  title  information,  up  to 

80  characters* 

2 

13 

K 

No.  of  Lamdaa  (wavelengths)  to 
be  provided  (a  maximum  of  185; 
see  card  10) . 

3 

F5.3 

EPS 

E  from  Program  OSSC. 

4 

12 

NPTS 

No.  of  08cillator8  (a  maximum 
of  8). 

5 

7F8.2 

FR(I) 

NU  values  from  Program  OSSC. 

6 

7F8.4 

S(l) 

S  values  from  Program  OSSC. 

7 

7F8.4 

G(I) 

GAM  values  from  Program  OSSC. 

8 

13 

M 

No.  of  AI's  (a  maximum  of  10) 
(see  card  11). 

9 

13 

T 

Temperature  of  the  object  which 
is  emitting  (in  degrees  K) . 

10 

10F8.3 

LAMDA(J) 

Lamda's  at  which  omittance 
polarisation  is  calculated  for 
each  AI(N). 

11 

F10.6 

AI(N) 

Angles  of  incidence  (1  per  card) 
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APPENDIX  D 

PROGRAM  LISTING  -  OSCILLATOR 

DIMENSION  RFC(  13, 600), 01(6025,02(602) 

DIMENSION  MINI  300) »RIN(  300) 

DIMENSION  LABEL! 6) 

DIMENSION  Lt5K(9),LBL(9),L8M(2) 

REAL  X S ( 25) ,H( 27) 

REAL  X(25),G(25),Q( 1000 )«A(20),DATA(250) 

OATA  LBL/36H  ALPHA  / 

OATA  LABEL/4H  R  , 4H  El  ,4H,c2  ,4H-£F  ,4H  N  ,4H  K  / 

C  . ••*•*•••*••« . . . . 

RE AC  I  5, 548 ) A, LIMIT, EPS, EST 
548  FORKAT(20A4,/I5,E15.5,F10.0) 

C  ♦*«••*•*...**• . . 

RE  AC ( 5 , 55 ) N 
55  FORMA) (15) 

C  . . . . . . 

REAC(5,5MXm,I*l,N) 

5  FORMAT (3F10.0) 

C  - 

WRITC(6,40)A,LIMIT,EPS,FST 

40  FOR/' A) (1H1,20A4,/'  LIMIT  OF  ITERATIONS  IS', 15,*  EPS  IS  ',£15.5,// 

ESTIMATED  FIX)*1  F12.5,/) 

C  - - - 

HR  I TE ( 6 , 470 ) 

470  F0RMAT(5X,' INITIAL  GUESSES  FOR  OSCILLATOR  PARAMETERS  ARE  •//  ) 
NM«N-1 

C  - - - 

WRITE(6,47l)  IXm,I«l,NM) 

471  FORMAT ( 14X , 'S' , 16X , 'NU* , 14X , 'GAMMA • , //( 5X, E12. 5,8X, E12. 5, 

.  6X.E12.5) ) 

C  - 

HR  IT E ( 6 , 472 )X ( N ) 

472  FORMAT  I 3X,  • £  I  INF  > -  •  ,F12.4) 

C  . . . . . . . 

RE AC ( 5 , 55 )  M 

C  - - - 

HRITE16, 1766IM 

1766  FORMAT ( '  CARO  COUNT  IS  »  ,15) 

C  - 

HR  ITT( 6, 2000) 

2000  F0RMATC1REFLECTIVITY  BASEO  ON  INITIAL  GUESSES'/) 

C  - 

HR  I  TE ( 6 , 487  ) 

C  . . 

RFaC(5,555HHIN(  I)  ,RINI  I ),  I»l,M) 

555  FORMAT (2F10.0) 

00  29  l»l,N 
29  X< I  )-SCRT(X(  I )  I 
2TI-0 

cc  io  j«i,r 

RI.N(J)»RIN(  J)/100. 

YY « f UNX ( N, X ,n I M  J  )  ) 
l-YY-RIN(J) 

2L  T»2L  T  +  2#2 

C  - 

HRlTE(Slt)  WIN(J),RIN(J),VY,Z 
10  COMInLl 

C  - - 

WKITJC,IU2)2LT 

1112  FORMAT (  5X, •  SUM  CF  THE  DIFFERENCES  SQUARED  IS  • ,612.4) 

CALL  FUNCTI I  M,  HI'/,  •<  IN ) 


| 

! 


V 

5 

•I 


« 


4 
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i  'X'-HiJS-T1'  -V; 


T.*-3«AS«i« 


»w*r»  rrr^aa-'-r^P’fVe’H^?  f’OTBI 


EXTERNAL  FUNCT 

IFILIMIT  . EU.  0)G0  TO  250 

CALL  FMFP! FUNCT, N,X,F,G,EST, EPS,LIMI T, IER*8) 

250  CONTINUE 

00  19  1*1,  N 

19  XSII)<XII)*t2 

C  - — - 

WRITE!  6,. 974 5) 

9745  FORMAT  1 1H1) 

C  - 

WRITE(6,473)F, IER 

.473  F0RMATI3X, ///‘FITTING  HAS  BEEN  COMPLETED  WITH  FIX)*  »»E15.8,// 

AND  ERROR  PARAMETER  -  ‘,I5//‘  FITTEO  VALUES  FOR  OSCILLATOR', 

.»  PARAMETERS  ARE  •/) 

C  - - - 

WRITE (6, 471) iXS(IJ, I*1,NM) 

WRITE(6,472)XSCN)  * 

WRITE! 6,487) 

487  FORMAT ! 1H1, 3X‘ ENERGY  RIEXP)  R(FITTEO)  DIFFERENCE1/  ) 

ZLT*0 

DO  13  J* 1 ,  M 
YY«FUNX<N,X,WIN! J)  ) 

Z* YY-R IN ! J ) . 

ZLT-ZLT*Z*Z 

C  - - - - 

WRITS(6,16)WIN!J),RINIJ),YY,Z 
16  FORMAT ( 1X.4F12.5) 

13  CONTINUE 


WRITE!6»1112)ZLT 

ENC 


SUBROUTINE  FUNCT I! M,W IN.R IN) 
DIMENSION  WIN! 300 ),RIN(300) 

RETURN 

ENTRY  FUNCnN,X,VAL,GRACJ 
REAL  X!25) ,GRA0!25) 

REAL  W(10),S(10),G(10) 

REAL  P ! 300  l  , Q ( 300 ) , V 1 300 ) 

COMPLEX  Z,0( 300,10),E! 300), Hi  300) 

NN*\/3 

DO  1  I-l.NN 

L«3* !  I-i) 

SI  I  )*X!L  +  1) 

W!  I  )  * X  ( I.  ♦  2  ) 

1  G! I  )*X ! L+3 ) 

E I  N  F  *  X  ( N  ) 

00  2  J*1 ,NN 
W2“l ,/W! J ) *»2 
C2-C! J)**2 

00  2  1*1,  M 

WW-WINI I  )*W2 

2  0! I,J)*SIJ) /CMPLX! l.-WW«*2.-WW*C2> 
00  J  I  *  I » M 

El  I )*E INF*»2 
00  4  J*  1 , N'l 
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4  £1 I ) =G ( I ) +S ( J ) *D( I t J ) 

E( » )*CEXPI . 5*CL0G( El  I)  ) ) 
IF(R£AL(E(I)).LT.C.)E(I)*-E(I) 

Pi I1*CABS(E< I 1 )«*2-2.*REALIECIl 1*1. 
‘  0( I i-P ( I 1+4.*REALIEI  1)) 

3  VI  I)~Pm/Q(I)-ftl*<I) 

VAL*0 

,00  33  1*1, M 
33  VAL*VAL*Vl I  1**2 


WR I TE ( 6 « 76 1  VAL 
76  FORMATI 1X.E14.5) 

00  7  1*1,  .M 

Z=CCNJG(  E<  I  )  1 

7  HI  I )=8.*V<  I )*.5*t Cl  1 1*1 Z-l. 1-PI  I l*IZ+i. J 1/El I l/QI I )**2 
00  6  J*1 ,NN 

L*  3* ( J-l 1 
GR AC  I L» 1 1 *0 
GRACIL+21-0 
GRACI L+3}*0 
, W2*G I J 1/WI J 1 **2 
W3*-C( J)**2/WI Jl**3 
W5*-2./WI Jl**5 
DO  6  I » 1 , M 

GRACIL*1)*GRACIL+ 11+REALiHl 11*01  I, J) 1 
Z-H  «  11-01  I , J 1 **2*V INI  I  1 

GRACIL*21aGRAO(L*21*REALIZ*GMPLXlWINII 1*W5,K31  1 

GR AC (L*3)*GRAD(L*3J— AIMAGIZ1 *W2 

CGNTINUE 

GR  AC  I N  1  *0 

CO  6  I  *  1 ,  K 

8  GRACI N 1 =GRAD<  N 1+REAL I  HI  1 1 1 *EINF 
RfclURN 

END 


FUNCTICN  FUNX I N, X , W I N 1 
REAL  XI25) 

REAL  WdOl.GIlOlfSIlO) 

NN-N/3 
DO  1  I  *  1 » NN 
L«  3* I I-l 1 
SI  I  )*X(L*1) 

Wl  Il-XIL+2) 

G(  n*XIL*3) 

CONTINUE 
COMPLEX  E 
E*XIN)**2 
DO  2  J*  l , NN 

E«E*SiJ>*«2/CMPLX< l.-IWIN/W J l**2)**2,-WIN*CI J1**2/WI J)**21 
CONTINUE 

E-CEXPI ,5*CL0GIEJ 1 
IFIREALIEJ.LT. 0. IE *-E 
FUNX-CA8SII E-1.)/IE*1.) )**2 
RETURN 
END 
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APPENDIX  E 

PROGRAM  LISTING  -  EMISPOL 

DIMENSION  LAMDA( 135) ,AM( 185),AX( 185) , A I( 10) , ER ( 185 )  ,EL ( 185  ) 
DIMENSION  0E( 185) , NR( 185 ) ,NLI 185) ,DNI 185) • I TI TLE 1 10) .NSKYI 9, 185) 
DIMENSION  FR(6),S(3),G(S),EE(185) 

DIMENSION  PERPCU195) 

COMMON  FR, S * G  " 

COMMON  NPTS.EPS 
REAL  LAMDA. NSKY,NTAR,NR,NL 
INTEGER  T 
INTEGER  SKY 
C  INPUT 

100  REAC  (5.5)  IT  I TLE 

5  FORMAT (10A4) 

WRITE  (6,6)  ITITLE 

6  FORMAT  (1H1.10A4) 

REAC  (5,10)  K 

10  FORMAT  ( 13) 

REAC  (5,12)  EPS 
12  FORMAT  (F5.3) 

REAC  (5,21)  NPTS 
21  FORMAT  (12) 

REAC  (5,20)  (FR(I), 1=1, NPTS) 

20  FORMAT  (8F8.2) 

REAC  (5,25)  (S(I), 1=1, NPTS) 

25  FORMAT  (9F3.4) 

REAC  (5,25)  <G(I), 1=1, NPTS) 

REAC  (5,40)  M 
40  FORMAT  (13) 

56  REAC  (5,60)  T 
60  FORMAT  (13) 

609  REAC  (5,610)  ( LAMCAI I ) , 1=1, K ) 

610  FORMAT  (10F3.0) 

N*0 

645  N* N ♦ l 

REAC  (5,650)  MINI 
650  FORMAT  (F10.6) 

A I ( N )  =  90 .0- A  I ( N ) 

DO  660  1*1, K 
NSKY ( N  ,  I ) =0 • 0 
IF  (N.GE.M)  GO  TO  600 
GO  TO  645 
660  CONTINUE 
C  CALCULATIONS 
600  00  700  JJ»l,K 

CALL  ANAK(LAMOA( JJ),AN( JJ),AK( JJ5 ) 

700  CONTINUE 

DO  200  1 1  =  l  ,  M 

AI  I  II)=AI( II )*3. 1416/180. 

00  150  JJ«i,K 

Q»SCRTC4.*UN(  JJ)**2)*(AK(  JJ)»*2)*(  ( AN  ( J  J )  **2  )  -  ( AK  (  J J ) **2 ) -S I N ( A I ( 
III  )  >*<*2)**2> 

A=Q.5*(+AN( JJ)**2-AK( JJ )**2-(SIN(AI( 1 1 ) ) )**2+Q) 

B=0.5*(-A.N(  JJ)**2*AK(  J  J  )  **2*  (  S I  N(  A I  (  1 1 ) ) )  **2+Q ) 

RR  =  ( ( SCRT( A l-C05( A I ( II)) )«*2  +  B )/  ” 

1  ((5CRT(A)+C0S(AI(II)))#*2+8) 

RL  =  (P.R«(  (SQRT(  A)-SIN(AI(  l!) > *TAN( A I ( 1 1 )) ) «*2*B )  J / 

1  ( ( SORT ( A)+SIN(AI(II))«TAN(AI(II) ) )**2*8) 

NTAR-i. l9E4/(( LAMCAI J J ) **5 ) * ( ( E XP ( l. 4388E4/ (LAMOA ( JJ) *T) ) )-l.O) ) 

ER( JJ)=l.O-KR 

ELI JJ)=l.0-RL 

EE ( JJ ) ».5<M  ER( JJ ) ♦ELI JJ ) ) 


Preceding  page  Mink 
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DEI JJI-ELt JJ)-ER| JJ) 

NR  I JJ)«|RR*NSKY( II ,JJ )/2.0>+(ERI JJ)*NTAR/2.0) 
NL(Jj)*(RL»NSKY(IItJJ)/2.0)+<EL(JJ)«NTAR/2.0) 

PNI  JJ)«NL(  JJ)-f<R(  JJ) 

PERPOL ( J J  ) ■ I  EL ( JJ ) -ER I J J ) ) / 1 EL I JJ ) *ER (JJ )  ) 

150  CONTINUE 

All II)«AI( I n # 180. /3. 1416 

WRITE  16,70)  AI(tI)fT,( LAMDAt  L) ,AN(L),AK(L) *NSKY( II,L),ER(L),ELIL) 
1  <0E  I L  t  ,  NRIL  ) ,  NL I L  )  ,  ONI  L  I , EEIL )  »  PERPOL  IL),Lal,K) 

70  FORMAT  l//»AI»  • , F3. 3/, • TEMP-  • , I  5,/ , IX , • LAMDA • ,3X , • AN* ,6X , • AK» ,6X 

1,  'NSKY*,7X,  '£R,,9X,,EL',9X,,0E,,9X»,NR,,9X,,NL'»9X,»DN,,9X»*£E**TX 

2, *PERPCL*/»  IF6.2,2X,F6.4,2X,F6.«,2X,E10.*,1X,E10.4,1X,E10.4,1X 

3,  E 10.4, IX, E9. 3, IX, E9. 3, IX, E9. 3, IX, El  0.4, IX,  E10. 4) I 
200  CONTINUE 

GO  TO  100 
END 


SUBROUTINE  ANAKILAMBOA.AN.AK) 

DIMENSION  FR(8),S(8),G(6) 

COMMON  FR.S.G 
COMMON  NPTS.EPS 
REAL  lambda 

F=lOOOC./LAf1RCA 

AE*C. 

Btt-O. 

DO-  110  KK-l.NPTS 

00*  (SIKK)«(  IFRIKK) )  *<*2 ))  / 1  II  FR(  XK  >>**2-F*»2 )  *»2*  C  G  IKK )  *FRI  KK)  »F)  ** 
12) 

AE«CQ* I IFRIKK) )**2-F**2)*AE 
8B«  I *  5 "G I KK ) *FR I KK ) • F • Q  C ) *80 
110  CONTINUE 
AA-AE+EPS 

AN-SQRTI .5*1 AA+SQRT I AA««2*A. •109**2)  )  )  ) 

ak=hb/an 

RETURN 

ENC 


APPENDIX  F 

PROGRAM  LISTING  AND  INSTRUCTIONS  FOR  IBM  SCIENTIFIC 
SUBROUTINE  -  FMFP 


(This  material  in  this  Appendix  is  a  direct  copy  cf  an 
IBM  Scientific  Subroutine  as  it  appears  in  IBM  Manual 
GH 20-0 20 5-4  of  System/360  Scientific  Subroutine  Package, 
Version  III,  Programmer's  Manual.) 
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Extremum  of  Functions 


Subroutines  FMFP  and  DFMFP 


These  subroutines  perform  the  calculation  of  an  un¬ 
constrained  minimum  of  a  function  of  several  vari¬ 
ables  using  a  method  proposed  by  Davidson.  The 
underlying  method  is  described  in  the  article  by 
R.  Fletcher  and  M.J.D.  Powell,  "A  Rapidly  Con¬ 
vergent  Descent  Method  for  Minimization",  Computer 
Journal,  vol.  6,  iss.  2,  1963,  pp.  163  -  166. 

It  is  assumed  that  the  function  f  of  the  n  variables 
xi,  ...  ,  xn  (abbreviated  as  argument  vector  x)  may 
be  computed  together  with  its  gradient  vector  g(x) 
for  any  point  x.  The  generalized  Taylor  expansion 
for  functions  of  several  variables  is 

1  T 

f(x+u)  *  f(x)  +  g(x)  •  u  +—  u  G(x)u  +  higher  terms 

where  g  is  the  gradient  vector  and  G  the  matrix  of 
second  order  partial  derivatives.  Vectors  are  as¬ 
sumed  to  be  column  vectors;  u^  means  transpose  of 
vector  u.  It  is  assumed  that  in  the  neighborhood  of 
the  required  minimum  xmjn  the  function  is  approxi¬ 
mated  closely  by  the  first  three  terms  of  its  Taylor 
expansion,  giving 


f(x)  *  f(x  ,  )  +-r  (x  -  x  )  TG(x  .  )(x-x  ) 

'  '  min  2  '  min  '  min  min' 

since  g(xmm)  *  0.  Then  the  gradient  is  seen  to  be 
approximately  g(x)  *  G(xmln)  (x  -  xrata). 

Assume  now  that  the  symmetric  matrix  G  is 
positive  definite.  Then  the  following  equation  holds 
true: 


min 


,<rl<W 


g(x) 


which  would  allow  xmjn  to  be  calculated  in  one  step 
if  G~l  (xmin)  were  available. 

To  approach  G'1  (xmin),  a  method  of  successive 
linear  searches  in  G-jon jugate  directions  is  used. 
Starting  with  the  identity  matrix  G^°)  *  I,  a  sequence 
of  symmetric  matrices  G^1)  is  generated  which  tends 
to  G“*.  At  the  (l+l)8t  iteration  atet>  a  linear  search 
is  made  in  direction  hi1)  «  -  Gl1)  g'1),  where  gl1)  is 
an  abbreviation  for  g(x(0).  By  means  of  the  linear 
search  the  minimum  of  y(t)  ■  f(x0)  ♦  t-M1))  is  deter¬ 
mined,  giving  argument  xi1*1)  "xl1)  +  tj  •  hi1). 

The  argument  of  the  minimum  x(1+1)  on  the  line 
through  xl‘)  in  direction  hi1)  is  determined  by  the 
relation  that  scalar  product  (gl1+1),  hi1))  ■  o. 


New: 


and: 


;<»>., 0). 

n-1 

E 

tt  h<‘> 

i-j 

0>. 

n-l 

E 

VGh<‘> 
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Therefore: 


coaler  product  (g(n),  h^)  -  £  t  <Oh(i>,  h0)) 

i-J+1 

Suppose  now  that  the  vectors  h(°),  hM,  ... , 
are  G-conjugate,  satisfying  (GhO),  hO))  « 

0  for  i  f  J.  Then  (g(n),  h(J))  =  0,  and  since  h(°), 
hM,  ...»  h(n-1)  form  a  basis,  g(n)  *>  0  and  x(n)  « 
xmin.  This  shows  that  the  minimum  is  located  at 
the  n*h  iteration  for  a  quadratic  function  when  using 
successive  linear  searches  for  G-conjugate  direc¬ 
tions. 

For  the  generation  of  G-conjugate  directions, 
start  with  h^  »  -g^  and  calculate  successive 
directions  h^  by  means  of  =  -G^g®\  where 
G^  is  modified  to  G^  +  ^  so  that  h^  is  an  eigen¬ 
vector  of  the  matrix  G^+  ^  G  with  eigenvalue  1.  This 

ensures  that  G^  approaches  G  1  as  approaches 

x  ,  .  An  easy  calculation  shows: 
nun 

(Ul)  (1)  dx  ■  dxT  G^dg  •  dgTG(l) 

dxT.dg“  dgV>dg 

with  dg  «=  g^  *  ^  -  g^ 
dx-x^-x^ 


where  all  vectors  are  regarded  as  column  vectors, 
and  superscript  T  means  transpose  of  column 
vector— that  is,  row  vector. 

The  strategy  adopted  for  termination  of  the  suc¬ 
cessive  linear  searches  is  as  follows: 

1.  If  the  function  valuo  has  not  decreased  in  the 
last  iteration  step,  the  search  for  the  minimum  is 
terminated  provided  the  gradient  is  already  suffi¬ 
ciently  small;  otherwise,  the  next  step  la  in  the 
direction  of  stoepost  descent. 

S,  If  the  argument  vector  and  the  direction  vec- 
changa  by  vory  small  amounts,  and  at  least  n 
Iterations  are  performed,  the  minimisation  is  ter¬ 
minated  again. 

3.  if  the  number  of  itorationa  exceeds  an  upper 
bound  furnished  by  ths  user,  further  calculation  is 
bypassed,  and  an  error  code  la  set  to  1  indicating 

poor  convergence.  Mathematics — Extremes  of 
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4.  If  one  of  the  successive  linear  searches  indt- 
cates  that  no  constrained  minimum  exists,  further 
calculation  Is  bypassed  again,  and  the  error  code  is 
set  to  2  Indicating  that  It  is  likely  that  no  minimum 
exists. 

The  i**1  term  G'lJ  is  reset  to  the  identity  matrix  if 
there  is  Indication  that  the  current  g(*)  is  not  posi¬ 
tive  definite,  or  if  the  formula  for  G@  +  breaks 
down  due  to  zero  divisors. 

The  linear  search  technique  used  in  subroutines 
FMFP  and  DFMFP  is  as  follows.  For  a  given  argu¬ 
ment  vector  x  and  a  vector  h  defining  a  direction 
through  x,  a  local  minimum  of  the  function  y(t)  » 
f(x+th)  is  determined.  This  means  that  a  value  ^ 
must  be  determined  such  that  y*  (tm)  =  0. 

Given  y* <t)  «  scalar  product  (g(x  +  th),  h),  there¬ 
fore  y(t)  and  y'(t)  may  be  calculated  for  any  value  of 
t,  and: 


y(0)  =  f(x)  and  y'(0)  *  (g(x),  h) 

In  case  y'(s2)  =  o,  tm  is  set  equal  to  82  and 
xm  m  x  +  s2,  h  is  used  as  argument  of  a  constrained 
minimum  on  the  line  through  x  with  direction  h.  In 
the  second  and  third  case  a  minimum  lies  between 
the  points  xl  =  x  +  sl-h  and  x2  «  x  +  s2  •  h/  that  is, 
tm  must  be  in  the  interval  (si,  s2). 

The  argument  of  the  minimum  is  obtained  by 
means  of  cubic  interpolation  using  the  function  and 
derivative  values  at  the  points  xl,  x2.  Let  X3  » 
x  +  s3  •  h  be  the  argument  of  the  minimum  of  the 
third  degree  interpolation  polynomial.  Then: 


■  s2  -o(s2  -  si)  >=sl  +  (a-  1)  (si  -  s2) 


with: 


y»  (s2)  +  w  -  z 
*  *  y  02)  -  y'(sl)  +  2w 


and:  a  *  3  ^jfaTs'l^  ♦  y'  (■*)  +  y’  (•*) 


w  «J[*8  -y' (■  i)  *  y*  («2)] 

If  f(x3)£»  f(xi)  and  f(x3)  $f(x2),  xm  is  set  equal  to 
x3  and  used  as  argument  of  the  wanted  minimum 
along  the  given  line.  Otherwise  the  interval  Qd,  x2) 
is  reduoed  by  replacing  xl  by  x3  if  f(x3)  $  f(xl)  and 
(g(xl),  h )  <  0,  and  by  replacing  x2  by  x3  la  all  other 

cases.  Then  the  interpolation  process  is  repeated  Mathematics — Extreaua  of 

for  this  new  reduced  interval.  Functions 
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"BURR OU  T"I  ^TThTP 

purpose  ‘  .  . . 

TO  FIND  A  LOCAL  MINIMUM  OF  A  FUNCTION  0E__3EVCRAI._VARIABLES 
BY  THE  METHOD  OF  FLETCHER  'AND' POhELL 

"usage  — 

C  A  L  LJLH  FH(  FU  N  C  T  fN.X  ,F  »  G.E  S  T  ,E  PS  ,  LIMIT, IER, H) _ 

DESCRIPTION  OF  PARAMETERS  _ _ _ _ 


FMFP  1« 
.FHFP  2d" 
FMFP  10 

7hfP  «5“ 
FMFP  50 
FMFP  6<F 
FHFP  70 
FMFP  60' 
_F*FP  90 
FHFP“l00 
FHFP_t  10_ 
FMFP  120 
FMFP  110 


FUNCT 


-USER-WRITTEN  SUBROUTINE  CONCERNING  THE  FUNC T I ON  fO  FHFP  lOO 

RE  MINIMIZED.  IT  MUST  RE  OF  THE  FORM _ fMFP_150 

SUPHfiuTINE  FUNCT (N*ABG , V *L , GRAD )  FHFP  160 

AND  MUST  SERVE  THE  FOLLOWING  PURPOSE  _  _FMFP_170. 

FOR  EACH  N-DIMEN3I0NAL  ARGUMENT  VECTOR  ARG,  FHFP  180 

FUNCTION  VALUE  ANO  GRADIENT  VECTOR  MUST  RE  COMPUTEPFMFP  190 

'  AND,  ON  RETURN,  STORED  IN  VAl  AND  GRAD  "RESPECTIVELVFMFP  200 

-  NUMHFH  UF_VARIABI,ES  _  _  _ _ FMFP_2lO_ 

VEC 1  OR- o'F  DIMENSION'  N "CONTAINING  the  INITIAL  FmFP  220 

ARGUMENT  MHERE  THE  ITERATION  STARTS.  ON  RETURN,  FHFP  230 

X  HOLDS  THE  ARGUMENT  CORRESPONDING  to  THE  FHFP  2N0 

COMPHTEO  MINIMUM  FUNCTION  VALUE  FMFP  250 

-SINGLE  VAHIaRLE  CONTAINING  THE  MINIMUM  FUNCTION  FHFP  260 
value  on  RETURN,  I.£.  FaF(X),  _  _  _ :_f  MFP_. 

-  VECTOR  UF  DIMENSION  N  CONTAINING  THE  GRADIElilf  F«FP  280 

VECTOR  CORRESPONDING  10  THE  MINIMUM  ON  RETURN,  FHFP  290 

I.E.  G*G(X>.  FHFP  300 

-  IS  AN  FSTIMATE  OF  THE  MINIMUM. FUNCTION  VALUE.  FMFP  310 

-  TESTVALUE  REPRESENTING  THE  EXPECTtP  ABSOLUTE  ERROR. FHFP  320 


_  _ A  REASONABLE  CHOJCF  IS  10**(-6)j_1.E.  __  _ _ 

•30**E HH A I  GREATER  than  10**  C-D) ,  NmErE  D  IS  THE 
NU’*BFR  t)F  SIGNIFICANT  DIGITS  IN  FLOATING  POINT 

representation. 

LIMIT  -  MAXIMUM  NumhER  OF  iterations,  . 

IER  ~  -  ERRUP  PARAMETER 

_ IE»  a  OCEANS  CONVERGENCE  MAS  OBTAINED 

lEP'a  1  means  NO  CONVERGENCE*  IN  LIMIT  ITERATIONS 
IER  a- 1  meAnS  ERRORS  IN  GRADIENT  CALCULATION 
IER  ■  2  MEANS  LINEAR  SEARCH  TECHNIQUE  INDICATES 
IT  IS  Ll«ELY  THAT  THERE  EXISTS  NO  MINingM, 

H  -  MORMNG  STORAGE  OF  DIMENSION  N*(N*7)/2, 


REMARKS 


_FMFP  330 
"FHFP  390' 
FHFP  350 
FHFP  3b0 
FHFP  370 
FHFP  300 
FHFP  390 
'FHFP '900' 
FHFP  «10 
'FHFP  920 
FHFP  9i0 
FHFP  990 
FMFP  9*0 
F HF 0  960 


T)  T“E  SUBROUTINE  NAME  REPLACING  THE  DUMMY  ARGUMENT  FUNCT  FMFP  970 


must  S£  DECLARED  *s  extfrnal  IN  THE  CALLING  program. 

IT)  I *R  IS  SET  To  2  IF  ,  STFPPING  IN  ONE  OF  !«£  COHPOTED 
DIRECTIONS,  THE  FUNCTION  WILL  NEVER  INCREASE  WITHIN 

A  tolerable  m*f.ge_of  ARGUMENT,  _ 

IFR  a  2  HAY  OCCUR  ALSO  IF  TH£  INTERVAL  WHERE  > 
INCREASES  IS  SHALL  AND  THE  INITIAL  ARGUMENT  HAS 
RELATIVELY  FAR  Am AY  FRQM  THE  NINJMUM  SUCH  THAT  THE 
MINIMUM  has  OVERLEAPED.  THIS  IS  DUE  TO  THE  SEARCH 
TECHNIQUE  NHICH  DOUPLES  THE  STEPftlZE  UNTIL  A  POINT 
is  EUUND  WHERE  THE  FUNCTION  INCREASES. _ 


SUOROUTINES  and  FUNCTION  SUBPROGRAMS  REQUIRED 
FUNCT 


FMFP  980 
FMFP  990 
FHFP  300 
FMFP  510 
FHFP  520' 

FHFP  53" 
FHFP  590 
FHFP  350 
FHFP  560 
FHFP  3T0 
FHFP  38 O' 
FHFP  390 
FHFP  690 


61 


M 

C 

FMFP  610 

62 

C 

METHflO 

FNFPf  620 

63 

C 

THE  METHOD  IS  0E8CMISE0  IN  THE  FOLLOMlNO  ARTICLE 

FMFP'  630 

64 

C 

K.  FLETCHER  AND  H.J.Q.  PORELL.  A  RAPI&  DESCENT  nEthITd  fOlP'" 

FMpfT  640 

65 

c 

HTNTMf 2ATIQN. 

PMF*  650 

66 

c 

COMPUTER  JOURNAL  VOL. 6,  183.  2,  1963.  PP. 163-168. 

FMFP  660 

67 

c 

FMFP  670 

68 

c 

.FMFPf  680 

69 

c 

FMFP  691. 

70 

SUBROUTINE  FNFP(FUNCT,N,X,F,6,£3T#EP3,LlMlT»lER,H> 

FNFP  700 

71 

c 

FMFP  710 

72 

c 

DIMENSIONED  DUNHV  variables 

FMFP  720 

73 

DIMENSION  H(l).Xd).GCl) 

FMFP  730 

78 

c 

FMFP  740 

75 

COMPUTE  FUNCTION  VALUE  AND  GRADIENT  VECTOR  FOR  INITIAL  ARGUHENTFMFP  750 

76 

C*LL  FUNCT(N,*,F,6) 

FMFP  760 

77.  . 

C 

FMFP  770 

78 

c 

RESET  ITERATION  COUNTER  AND  GENERATE  IDENTI TV  MATRIX 

FMFP  780 

79 

IER»0 

FMFP  790 

«0 

KDUNT*0 

FMpp  800 

61 

NPsNtN 

FMFP  610 

62 

N3*N2*N 

FMFP  820 

63 

N3l*N3f1 

FMFP  830 

64 

t  K«N31 

FMFP  840 

65 

DD  4  J»t.N 

FMFP  850 

66 

H(K )*1 • 

FMFP  660 

67 

NJcN-.T 

FMFP  670 

68 

1F(NJ)5,5.2 

FMFP  680 

69 

2  DD  3  L«i.NJ 

FMFP  690 

90 

KLsKtL 

FMFP  900 

91 

3  H(KL)»0.  _  . 

FMFP  910 

«* 

4  K*KLt  1 

FMFP  920 

93 

_c 

• 

FMfP  939 

9« 

c 

START  ITERATION  loop 

FMFP  940 

95 

5  KOuNTbKDUNI  ♦  ! 

FMFP  950 

96 

c 

FMFP  960 

97 

c 

SAVE  FUNCTION  VALUFj  AgGUHtNT  VICTOR  AND  GRADIENT  VECTOR 

FMFP  970 

98 

oldf»f 

FMFP  980 

.99 

00  9  Jai.N 

FMFP  990 

too 

K*N*J 

FMF  P 1000 

t«l  ... 

H(K)*6(I> 

FMFP1010 

102 

K*K+N 

FMFP1020 

1«3  . 

M(K)*X( t) 

FMFP1030 

104 

c 

F  MF  P 1 040 

ins _ 

_CL 

DETERMINE  DIRECTION  VECTOR  H 

FMFP1050 

106 

49JAN3 

FMFP1060 

107 

1*0. 

F  MF  PI  070 

108 

00  8  L«1»N 

FMFP1080 

109 

T*T-GfL)*H<«) 

F  MfPl 090 

no 

IP(L-J)6»7,7 

FMFP1100 

111 

6  K*K  *N— L 

fmfpiuo 

112 

GO  TO  8 

FMFP1120 

113  . 

7  K«K*l 

FMFP1130 

116 

6  CONTINUE 

FMFP1140 

US 

9  H( J)*T 

FMFP1 150 

116 

c 

FMFPU60 

117 

£_ 

JMFP1U0 

118 

OT*0. 

FMFPltlO 

119 

HNRMaO. 

FMFP1I90 

ito  '  gnrm«o.  ""  fmfpizoo 
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1?1  C  _  _  FMFP1210 

122  C  CALCULATE  DIRECTIONAL  DERIVATIVE  AND  TtSTVALUES  FOR~ DIRECTION  FMFRI22D 

123  _ C _  VECTOR  H  ARP  GRADIENT  VECTOR  6. _ FMFP1230 

120  _  ~00  10  J»1,N  "  FMFP12O0 

125  HNRM«mnRMaASSIM(J1)_  _  _  _  _  _  _  FMFP12S0 

126  GNRMsGNRMAASS(G(J>)  ~  ”  '  '  F*FP1260 

127  10  OY*PY*MfJ)oGU)__ _ _  __  _ _ _  FMFP1270 

128  C  '  FMFP12B0 

129  _  C _ _ REPEAT  3FAPCH  IN  DIRECTION  OF  STEEPE.ST  DESCENT  IF  DIRECTIONAL  FHFP1290 

130  C  DERIVATIVE  APPEARS  TO  RF.  POSITIVE  OR  ZERO,  FMFP1300 

m _ _  JF(0Y)n#5|.51  _ _ _ _ FMFPU10 

132  C'  "  FMFP1320 

133  C  _  REPEAT  SEARCH  IN  OIRECTlONjOE  STEEPEST  DESCENT_JF  DIRECJION _ FMFP1330 

134  C  VECTOR  H  IS  SHALL  COMPARED  TO  GRADIENT  VECTOR  6.  FMFP13R0 

1 35  _ 1 1  1E(HNRM/GNRM~FPS )51.5 1.12 _ _ _ F MFP 1 350 

136  C  FMFP1360 

137  C  SEARCH  MINIMUM  ALONG  DIRECTION  H_ _ FHFP1370 

138  C  .  FHFP1380 

139  C  SEARCH  ALONG  H  FOR  POSITIVE  DIRECTIONAL  OERJVATIVE _ FHFP1390 

100  12  FY»F  FPFP1O00 

.101 _ ALFA«?.»tEST-F)/OV _ FMFP10JO. 

102  AM0DA> 1 .  FMFP1020 

10J  _  C  FMFP1O30 

104  ~  C  USE  ESTIMATE  FOR  STEPS IZE  ONLY  IF  if  IS  POSITIVE  AND  LESS  THAN  FMFP1000 

105  C  1,  OTHERWISE  TAKE  1.  AS  STEPSIZE _ FMFPlOgO 

106  IF(ALFA) 15,15,13  FMFP1060 

107  _ 13  IF(AlFA«Am80A110»15«1R _ ?HFPl«70 

100  14  'amboa*aLfa  "  FHFP1O0O 


15  ALFA«0. 


SAVE  FUNCTION  AND  DERIVATIVE  VALUES  FOR  OLD  ARGUMENT 
16  FX»FY 

_ OVsDY  _  _____ _ 


STEP  ARGHMFNT  ALONG_H  _ 

DO  17  I«1,N 

17  X(I)»X(T)*A«0nA<*HCI) _ 


_  CQMPIITE’FUNCTION  value  AND  6PADIENT  FOR  NEW  ARGUMENT _ 

CALL  FUNCT(N,X,OG) 

..  fv«p  .  . _ _j _ _  _ _ _ _ 

CCJMPIITF  DIRECTIONAL  DERIVATIVE  DY  FOR  NEN  ARGUMENT,  TERMINATE 

SEARCH,  TP  DY  IS  POSITIVE,  IF  DY  IS  ZERO  THE  MINIMUM  IS  FOUND 

.  DY«0,  _ 

DO  18  I«l,N 

is  Dv«DYfG(n*Mm 
IP(Dy) 10,36,22 

term t na if  search  also  if  the  function  value  indicates  that 

A  MINIMUM  HAS  SEEN  PASSED 

19  IF(FY-FV)20, 22/22 

REPEAT  search  AND  DOUBLE  STEPSIZE  for  FURTHtR  SEARCHES 

20  AM0DAOAMBDAAALFA 
OLFAaAMKOA 

End  OF  SEARCH  LOOP _ _ , _ 

terminate  IF  THE  CHANGE  In  ARGUMENT  GETS  VERY  LAIIGe 
lE(HN»M*AMBDA-l.E10)l6,t6,2t 


FMFP109D 
FMFP1500 
FMFP1510 
FMFP1520 
FMFP1530 
>MFPl50b 
FMFPI550 
FMFP1560 
FMFP1570 
FMFP15S0 
FMFP1590 
F  MFP1 600 
FMFP1610 
FMFP1620 
FMFP1630 
FMFP1640 
FMFP1650 
FMFP1460 
FMFP1670 
FME«M6S0 
FMFP1690 
FMF  P 1  TOO 
FMf  *1710 
'FMFP1720 
FMFM7AD 
FMFP1740 
FMFP1750 
EMEP1760 
EMEP1770 
FMFP17S0 
FMFP1790 
EMFP1S00 
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Mi 

C 

3 

F  MM»  1 0 1  ft 

Mi 

c 

liniik  of arch  technique  indicates  that  no  minimum  exist# 

FMFP1B20 

Mi 

21 

If  R«2 

FHFPM30 

M4 

retuhw 

FMFPI84<r 

MS 

c 

FMFP1«5« 

M6 

c 

INTERPOLATE  CUBICA'I.LV ‘IN  THE  INTERVAL  DEFINED  8V  THE  SFARCH 

FHFP1060 

M7 

c 

AbOVF  A4ft  COMPUTF  THf  ARGUMENT  X  FOR  WHICH  THE  INTERPOLATION 

FMFP1B70 

MB 

c 

polynomial  is  minimized 

FMFP1060 

MB 

22 

T*0. 

FNFP1BB0 

MO 

23 

IF(AHB0A)?<i,36,?4 

FMFPMOO 

Ml 

24 

Z*3.*(FV-FV)/AHRDAtftXADV 

FMFPM10 

M2 

ALF AsAHAxi ( ABS(Z) * ARS(DX) # ABS(DV) ) 

FMFP1B20 

Mi 

d*lfa»z/alfa 

FHFPM30 

1B4 

DAlFA*DAlF A*DALFA-DX/ALFA*0V/ALFA 

FMEP1B40 

M5 

IF(PAI  FAJSl. 25,25 

FHFP1B50 

Mfc 

25 

W«ALFA*BORT(DALFA) 

FHFP1B60 

M7 

*LF  AsDy-OVtW+W 

MO 

IF  (ALFA)  250*251 *250 

MB 

259 

ALF A*(0V«7fW) /ALFA 

200 

GO  TO  252 

201 

251 

ALFA*(/*OV«W)/(Z*DX*ZtDV) 

202 

252 

alfasalfawamoda 

203 

DO  26  1*1*4 

FMFPlBOft 

204 

26 

X(I)3X(m<T-ALFA)*H(I) 

fmfpmbo 

205 

c 

FMFP2000 

206 

c 

TERMINATE*  m  the  value  of  the  ACTUAL  FUNCTION  AT  X  is  LESS 

FMFP20M 

207 

c 

than  the  FUNCTION  values at  the  INTERVAL  ENDS.  OTHERWISE  REDUCEFMFP2020 

200 

l 

THF  Interval  BV  CHOOSING  one  END-POINT  equal  TO  X  AND  REPEAT 

FHFP2030 

200 

c 

thf  interpolation,  which  end-point  to  choosen  depend#  on  the 

FMFP2040 

210 

c 

value  of  the  function  and  its  gradient  at  X 

PMFP2050 

211 

c 

FHFP2ft60 

212 

CALL  FUNCT(N,X,F,G) 

FMFP2070 

213 

IF(F-FX12T*>7,20 

FMFP2000 

214 

27 

IF(F-FY)J4, 36,20 

FMFP20B0 

215 

20 

dalfa*o. 

FMFP2100 

216 

00  ?V  1*1,4 

FMFP2110 

217 

29 

0*LFA»0ALFA*Gm*M<U 

FMFP2l2ft 

210 

IF(OALFA)30, 33*33 

FMFP2130 

21B ... 

30 

IF(F-FX) 32, 31*33 

.FMFP21 40 

2?0 

31 

IF(O*-0*LFA)3P,36,32 

FMFP2I50 

221 

32 

F*«F 

FMFP2160 

222 

OVlDAl.F  A 

FMFP2170 

22i 

t*ALF  *  .  _  _ _ _ _ _  _ _ _ _  .  .  .  _ 

FMFP2100. 

224 

am«oa*alfa 

FMFR21Bft 

225.. 

GO  TO  25 

FMf P220D. 

226 

33 

1F(FY-P135,3'»,J5 

FMFP2210 

227 

34 

IF (OY-t) ALFA)  35,  36, 35  . .  . ... . . ..  ...  .  .  .  ....  _.  .  _ 

FMFP2220. 

220 

35 

FV«F 

FMFP223ft 

220 

DVaDALFA 

FMFP2240 

230 

amwoawambda-alfa 

FMFP2250 

231 _ 

GO  TU  2? 

.FMFP2765L 

2}2 

c 

FMFP2270 

2*3 

c 

TFRMJNATE,  IF  FU4CTIUN  has  Nyi.  DECREASED  DURING  LA#T.  ITERATION _ 

_ _  _  _ 

234 

36 

M  (OI.UF-FtEPS)  51*38*30 

235 

c 

.  _ _ 

236 

c 

J 37  . 

_  X 

TWO  consecutive  iteration# 

230 

230  . 

SO 

00.37  J«1»N 

X»NTJ 

240 

HfK)wG(J)-H(K) 
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m 

242 

243  _ 

24« 

2*5 

246 

2*7 _ 

2*6 

249  _ 

250 

251 

252 

253  . 

254 

255  _ 

256 

257  . 

258 

259 

260 

261 _ 

262 

263 

264 

265 

266 

267  _ 

266" 

269 

270 

271 

272  "" 

273 

274 

275 
27b 

277 

278 

279 
2*0~ 

2*1 

2*2 

2*3 

2*4 

2*5 _ 

2*6 

2*7 

2*8 

2*9 

2*0 
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_C 

C 

c 

c 


tfSM  +  K 


39 


40 

41 

_42 

43 


TEST  LENGTH  0*  ARGUMENT  DIFFERENCE  VECTOR  AND  DIRECTION  VECTOR 
IF  AT  LEAST  N  ITFRAT10N8  NAVE  *EEN  EXECUTED,  TERMINATE#  TF 

doth  are  less  than  eps 

1FR*0  _ _ _  _ 

IF(KQUNT-N)42,39,39 

T#0.  __  _ 

z*o. 

DO  40  J«1 ,  N _ _ _ _  _  _ 

MNTJ 

H«H(K)  _ 

K«KTN 

_T«T6APSfH(i<)J _ 

2«Z+N*HfK) 

IF(NN»H-EPSI4!,41,42 _  __  _  _ 

IF(T-Ff»«j)56,56,«2 

TERMTNA TF rT^  NUMBFR- OF“iTERA Tf ON8  MOULD  EXCEED  "LIMIT 
1F(KQUNT«LIHIT)43>50#50 _ 

PREPARE  IJPDATINOJJF  MATRIX  H_  _ 

ALF*«0. 

DO  47  J«1,N _ _ _ _ _ 

K*J*N3 

b«0.  _  _ 


44 

*5 

46 


47 


46 


49 


on  46  c*i , v 

XL  *NtL  _  _ _ 

n*m*H(KL)MOO 

IF(I.-J)44,45,45  _  _  _ 

*«KSN-L 

GO  TO  44  _ 

K»K  ♦  1 

CONTINUE  _ _ _ _ _ ___ 

K»N*  J 

ALFA«ALFA*4*N(K)  _ 

M(J)«W 

REPEAT  3 F A P c H  " In“DI RCCT iWij F~$ T FF^iT'DFS CE" N  f  \f  RESULTS' 

ARE  NQT  satisfactory 

IF(Z*ALFA14*, t,48 

UPDATE  MATRIX  M  ’ 

X«n31  _ _ _ _ _ 

DO  «9  L • 1 1  V 
KL*N2*L 
DO  49  J«L#  V 
N J«N2a J 

MfK)«H(K)*rtfKL)*HeNJ}/Z-H(L)*HiJ)/ALFA 
KBK  6  1 


FHFP2360 

FMFP2390 

FHFP2400 

FHFP2410 

FMFP2430 

F MFP2O40 

MFP2450" 

FMFP2460 

FHEP2470 

FMFP2480 

FHFP2490" 

FHFP2500_ 

FHFP2510 

FMFP2520 

FMFP2530 

FHFP2540 

FMFP2550 

FMFP2560 

"FRFP2570- 

FHFP2560 

FHFP2590 

FMFP2600 

FHFP2610 

FHFP2620 

FHFP2630" 

FHFP2640 

FMFP2660 

FMFP2460 

FHFP2670 

FMFP2660 

FMFP2690' 

FHFP2700 

FMFP2710 

FHFP27 20 

FHFP27iO 

FHFP2740 

FHFP2750 

FNFP2760 

FNFP2770 

FMFP2760 

FHFP27V0 

FMFP2*00 

FMf P2*l« " 

FMFP2*20 

FNFP2*3 0 

FMF  (>2*40 

FNFP2850 

FMFP2460 


2*2 

GO  TO  5 

FMFP2"70 

2*3 

C 

END  OF  ITERATION  LOO* 

FMFP2*80 

2*4 

C 

FMF  P2*90 

2*5 

c 

NO  CONVERGENCE  after 

UNIT  ITERATIONS 

FMFP2900 

2*6 

50 

ifr«i 

FMFP2*10 

297 

RFTURN 

FMFP2928 

296 

c 

FMFP2910 

249 

c 

RESTORE  OLD  VALUES  OF 

function  and  arguments 

FMFP294* 

3*0 

51 

DO  52  Jp!,N 

FMFP295* 
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t 
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301 

302 

303 

K«n2*J 

52  xrj)«M(«)  ""  '  . 

CALL  FUNCT(N,X,F,G) 

FHFP2«60 

FNFP2970 

FMFP2060 

304 

305 

306 

307 
306 
309 

C 

C 

c 

c 

_ £_ 

REPEAT  SEARCH  In  DIRECTION  OF  STEEPEST  DESCENT  IF  DERIVATIVE 
FAILS  TO  a*  SUFFICIENT. Y  SNAIL 

IF(GNRH-£PS155.5S.S3 

~  FHFP2990 
FNFP3000 
~ FMFPJ010 
FHFP3020 

TEST  FIR  REPEATED  FAILURE  OF  ITERATION 

FNFP3030 

FHFP3040 

310 

S3  IF(IER)56#54,54 

FMFP3050 

311 

54  IER«-1 

FMFP3060 

312 

GOTO  1 

FMFP3070 

313 

55  IFr«0 

FMFP3080 

31 9 

56  RETURN 

FNFP3090 

315 _ 

__ENp  ... 

X?lFP31ft5 

END  OF  FILF 


\ 


66 


